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ABSTRACT 

Magnetospheres of many astrophysical objects can be accurately described by the low-inertia (or 
"force- free" ) limit of MHD. We present a new numerical method for solution of equations of force- free 
relativistic MHD based on the finite-difference time-domain (FDTD) approach with a prescription for 
handling spontaneous formation of current sheets. We use this method to study the time-dependent 
evolution of pulsar magnetospheres in both aligned and oblique magnetic geometries. For the aligned 
rotator we confirm the general properties of the time-independent solution of Contopoulos et al. 
(1999). For the oblique rotator we present the 3D structure of the magnetosphere and compute, for 
the first time, the spindown power of pulsars as a function of inclination of the magnetic axis. We find 
the pulsar spindown luminosity to be L « (^ 2 fi^/c 3 )(l + sin 2 a) for a star with the dipole moment p, 
rotation frequency f2*, and magnetic inclination angle a. We also discuss the effects of current sheet 
resistivity and reconncction on the structure and evolution of the magnetosphere. 

Subject headings: pulsars - neutron stars - magnetospheres 
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1. INTRODUCTION 

Many astrophysical objects, including neutron stars 
and accretion disks, form highly magnetized magneto- 
spheres. Modeling of the structure of such magneto- 
spheres requires solving for the self-consistent behavior of 
plasma in strong fields, where field energy can dominate 
the energy in the plasma. This is difficult to do with the 
standard numerical methods for MHD which are forced 
to evolve plasma inertial terms even when they are small 
compared to the field terms. In these cases it is possible 
to reformulate the problem and instead of solving for the 
plasma dynamics in strong fields, solve for the dynam- 
ics of fields in the presence of conducting plasma. This 
is the approach of force- free relativistic MHD, or force- 
free electrodynamics (FFE; see, e.g., Komissarov 2002). 
In this Letter we present a new numerical method for 
solving the equations of FFE and apply it to the prob- 
lem of structure of pulsar magnetospheres. We calculate 
the shape of the magnetosphere for arbitrary magnetic 
inclination and determine the pulsar spindown power. 

2. NUMERICAL METHOD 

The equations of FFE can be derived as the low-inertia 
limit of MHD (Komissarov 2002). When plasma inertia 
and temperature are small, the balance of forces on the 
plasma is pE + ~ j x B = ( "force-free" condition) , where 
p and j are charge and current densities. When plasma is 
perfectly conducting and abundant to short out the accel- 
erating electric fields (E • B = 0), the Maxwell equations 
in special relativity together with the force- free condition 
give (Gruzinov 1999; Blandford 2002): 
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Eq. J2J) is a prescription for the plasma current that is 
needed to satisfy the constraints. The two terms have 
simple meaning: the first is the current perpendicular 
to fieldlines, given by the advection of charge density 
with the E x B drift velocity, while the second term is 
the parallel component of the current which reduces to 
(V x B) || in steady state. The system (llll't is hyperbolic, 
can be evolved in time, and is much simpler than the full 
MHD equations. 

In order to numerically solve this system we uti- 
lize its close relation to Maxwell equations and use 
the finite-difference time-domain (FDTD) method com- 
monly used in electrical engineering (Yee 1966). Its 
strengths are very low numerical dissipation and conser- 
vation of divergence-free condition to machine accuracy. 
Electric and magnetic fields are decentered component- 
wise on an orthogonal grid with electric components de- 
fined in the middle of cell edges and magnetic compo- 
nents in the middle of cell faces, allowing to compute the 
change of flux through cell face from circulation around 
cell edges. We use third-order Runge-Kutta time inte- 
gration, maintaining time alignment of the fields. To 
calculate the current in @ we use linear interpolation of 
the fields to the same locations on the grid (at the E- 
points). The second term in (JSJ) comes from the perfect 
conductivity constraint and is cumbersome to calculate 
numerically as it requires interpolation of both the fields 
and field derivatives. Instead of using the full eq. J5J) 
we solve an even simpler system: we advance Maxwell 
equations with only the perpendicular component of the 
current, and then after every timestep we subtract the 
accumulated E» component from the electric field to en- 
force E • B = 0. This is equivalent to replacing the paral- 
lel current with a resistive term <J\\E\\ (Komissarov 2004), 
where the parallel conductivity <7|| is taken to infinity. 

There are two conditions for the validity of the 
force-free approximation: electromagnetic energy den- 
sity vastly exceeds plasma pressure and inertia, and the 
plasma drift velocity is subluminal, so that E 2 < B 2 . 
The system lllll'l) does not enforce these conditions, and 
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initial force-free configuration may develop regions that 
violate one or both conditions. When E exceeds B, as 
may happen during the interaction of strong waves, the 
Alfven speed becomes imaginary and leads to an expo- 
nential instability, which physically would be averted by 
dissipation. Another example is the formation of cur- 
rent sheets, which is a generic feature of force-free evo- 
lution. In current sheets magnetic field can reverse and 
go through zero. This can violate both conditions, which 
means that the assumption of magnetic dominance is in- 
valid and one has to restore plasma effects such as pres- 
sure to sustain the current sheet. However, in the frame- 
work of FFE these effects are not included, so we try to 
keep current sheets unresolved on the grid, with a jump 
in fields over one cell. The FDTD method allows to sup- 
port such sharp discontinuities without spreading them. 
In order to have the method capture current sheets, we 
enforce the condition E < B after every timestep by re- 
setting the electric field to E = B in the regions where 
it exceeds B. The physics behind this drastic measure is 
related to dissipation processes that will happen if E x B 
drift approximation is violated and particles accelerate 
while decoupling from magnetic field lines. This dissipa- 
tion will restore E = B on plasma timescales; however, 
this physics remains to be clarified. We parametrize this 
physics with prescriptions for the dissipative current j± 
that is added to (0) when E > B. The prescription we 
use can be formulated as j± = a±E when E > B and 
jx = otherwise, where the perpendicular conductivity 
<j± — > co. Other prescriptions are also possible (Komis- 
sarov 2004). 

Our method is not the only way to solve FFE. By for- 
mulation our method is similar to Komissarov's (2004, 
2005), in that we also solve the non-conservative form 
of the equations with explicit resistive current prescrip- 
tions. The limiter procedure for enforcing E < B is simi- 
lar to what was used in McKinney's conservative scheme 
(2006a). The main difference, however, is that our FDTD 
discretization has inherently much lower numerical diffu- 
sivity. This is why our method is able to develop and 
sustain sharp current sheets unlike Komissarov's (2005) 
scheme, without the need to explicitly null the inflow ve- 
locity into current sheets as in McKinney (2006a). In 
fact, the location of the sheet does not have to be known 
in advance, which makes the method useful for evolution 
of complicated initial conditions. Below we describe the 
applications of our code to pulsars. 

3. ALIGNED PULSAR MAGNETOSPHERE 

Particle extraction from the stellar surface due to 
unipolar induction and subsequent pair-creation is likely 
to populate the pulsar magnetosphere with abundant 
nearly force-free plasma (Goldreich & Julian 1969). The 
detailed structure of the magnetosphere had remained a 
puzzle for a long time due to the difficulty of finding an 
analytical solution. The first numerical solution to the 
time-independent set of force-free equations (the "pulsar 
equation") was found by Contopoulos et al. (1999; here- 
after, CKF), and later improved by Gruzinov (2005a). 
The CKF solution consists of a region of closed field- 
lines extending to the light cylinder ( "closed zone" ) , and 
an "open zone" with asymptotically monopolar poloidal 
ficldlines. The requirement that the closed zone extend 
to the light cylinder was imposed from the outset, and 



was later relaxed by Goodwin et al. (2004), Contopou- 
los (2005) and Timokhin (2005), who find steady-state 
numerical solutions with closed zones of arbitrary ex- 
tent within the light cylinder. This raises the question 
which solution would be chosen by a time-dependent 
evolution. Recently, Komissarov (2005) used a time- 
dependent force-free code and a relativistic MHD code 
for this problem. While his force-free code had converged 
to a solution with all closed fieldlines, presumably due to 
excessive numerical diffusion, the MHD code produced a 
solution remarkably similar to CKF. The force- free code 
of McKinney (2006b) also converged to the CKF result. 
Generally, we are in agreement with these results, how- 
ever, intrinsic low dissipation nature of our code allows 
us to find features that were missing in prior studies. 

We discretize FFE equations in spherical coordinates 
on a uniform 1400 x 181 r — 9 grid with an absorbing zone 
to prevent reflections from the edge of the domain. The 
star has the radius i?* — 0.2Rlc, and the light cylinder 
Rlc = c/^* is at 300 cells. The magnetic field is ini- 
tialized to be a dipole B r = 2 fi cos 6 /r 3 , Bg = /is'mO/r 3 , 
Bcj, = 0, where /i is the magnetic moment. The electric 
field on the star is set to E = — SI* x r x B/c to sim- 
ulate a rotating conducting sphere, and the rotation is 
turned on as a step function at t — 0. During the initial 
evolution a torsional Alfven wave is emitted from the sur- 
face (Spitkovsky 2004, 2005). The wave transports space 
charge, and the associated electric fields set fieldlines in 
rotation. Near the star cancellation of waves from two 
hemispheres causes a growing region with zero poloidal 
current on closed fieldlines. Further from the star the 
fieldlines get stretched, and a pulse propagates along the 
equator. In ideal FFE the fieldlines cannot break or re- 
connect and the solution tries to stretch the fieldlines 
so they close at infinity. In the process the magnetic 
field strength on the equator drops and becomes smaller 
than the electric field, necessitating non-ideal physics, 
and eventually hits zero. This happens around 0.6i?£c 
after 1/4 of stellar rotation. Beyond this point the solu- 
tion forms a current sheet, supported by our resistivity 
prescription that keeps E < B. After 1 rotational period 
the magnetosphere forms the closed-open configuration, 
but with the Y-point at Q.GRlc- This is not the CKF so- 
lution, and in fact the configuration proceeds to expand 
the closed zone to Rlc by reconnection of the recently 
opened fieldlines. The rate of this expansion depends on 
the resistivity prescription and is 20 rotation periods of 
the star for our lowest resistivity model, but can be ac- 
celerated to 5 rotation periods if we allow for E to exceed 
B in the current sheet by limiting the resistive current. 

Once the closed zone reaches Rlc, the solution be- 
comes quasi-steady (Fig. la), and similar to CKF and 
McKinney (2006b). The fieldlines that reconnect outside 
Rlc (Fig. lb) break through emission of small plas- 
moids every 5 periods, dependent on resistivity, and the 
Y-point oscillates around (0.95 ± 0.05)Rlc- There is a 
sharp rise in the magnitude of the magnetic field on the 
equator at the Y-point (Gruzinov 2005a), followed by a 
jump by a factor of 10 2 to nearly in the current sheet 
(fig. lc, corresponding jump is a factor of 10 in fig. 3 
of McKinney (2006b), suggesting larger diffusivity) . The 
evolution towards the light cylinder after the initial for- 
mation of the current sheet takes much longer than the 
light cylinder crossing time, and the magnetosphere goes 
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Fig. 1. — Aligned dipole magnetosphere: a) Poloidal ficldlines of steady state solution. Thick line is the fieldline that touches the light 
cylinder, which is marked by a vertical line. Color represents toroidal magnetic field, normalized by [j,/(Rlc R*) 3 i b) Zoom in near the 
Y-point. Closed lines beyond the light cylinder are periodically thrown open by plasmoid emission; c) Magnitude of the magnetic field on 
the equator (in units of p,/R^ c ) during the evolution of the closed zone to the light cylinder. 



through a sequence of quasi-static equilibria similar to 
the states in Timokhin (2005) (see fig. lc). This evolu- 
tion will be missed if code diffusivity is increased. The 
spindown energy loss decreases as the closed zone ex- 
pands and saturates at (l±0.05)^ 2 f^/c 3 , consistent with 
previous works. 

4. OBLIQUE PULSAR MAGNETOSPHERE 

Real pulsars necessarily have nonzero angle between 
the magnetic and rotational axes (otherwise they would 
not pulse!), so understanding the structure of oblique 
pulsar magnetosphere is most interesting. Due to the 
lack of axial symmetry this problem has to be solved in 
full 3D. Encouraged by the convergence of our method 
to known solutions for the aligned rotator, we extended 
our code to three dimensions. Below we discuss gross 
properties of the field structure and energy loss. 

We use a Cartesian grid which avoids coordinate sin- 
gularities, but presents difficulties in imposing bound- 
ary conditions on a central sphere. Rather than us- 
ing stairstepping we force the fields to known val- 
ues inside the star with a smoothing kernel (similar 
to Komissarov 2005). Inside the sphere the field is 
set to B (t) = (3?/x(t) • f - f2(t))/r 3 , where = 
/i(sin a cos f2*i, sin a sin cos a), and a is magnetic in- 
clination angle. Electric field is forced to corotation val- 
ues. We use a grid of 700 3 cells with R* = 15 cells and 
Rlc — 40 cells. Such unrealistic ratio R*/Rlc is needed 
to resolve both the star and the wavelength 2itRlc on 
a limited grid. Without nonreflecting boundary condi- 
tions we only evolved the magnetosphere for 2 turns of 
the star, but this seems to be enough to establish steady 
state solution in the corotating frame. 

Fig. 2 demonstrates the force-free magnetosphere of 
a = 60° rotator after 1.2 turns of the star. Although 
the magnetospheric structure is intrinsically 3D, some 
insight can be gained from considering the shape of the 
ficldlines in the corotating frame in the plane defined 
by fi and 1~2» vectors. In this plane the fieldlines are 



reminiscent of the CKF solution with a closed and open 
zone and a current sheet. The main difference is that the 
current sheet oscillates about the rotational equator in a 
wedge with the opening angle of 2a and the wavelength 
of 1-kRlc- The ficldlines in this plane become straight 
beyond the light cylinder, corresponding to the inclined 
split-monopole solution found by Bogovalov (1999). For 
comparison, the fieldlines of a rotating dipole in vacuum 
look dipolar at all radii in /x-f2* plane, with no open 
zone or current sheet. In force-free solution the current 
sheet starts at the intersection of the light cylinder with 
the closed zone (fig. 2b) even in the oblique case. In 
the perpendicular plane the field is increasingly toroidal, 
reversing the sign in the current sheet (fig. 3a). 

An important question is the electromagnetic luminos- 
ity of the inclined rotator, since currently the vacuum 
formula L vac — 2/3^ 2 Sl*/c 3 sin 2 a is commonly used to 
infer the magnetic field of pulsars. Our 3D solution al- 
lows a direct measurement of this quantity. We measure 
the Poynting flux integrated over spheres of different ra- 
dius. After about 1 turn the solution settles to a con- 
stant energy flux, which we plot in fig. 3b as a function 
of magnetic inclination. We find that the formula 

-^pulsar = h — ^(1 + k 2 sin 2 a), (3) 

with coefficients k-y = 1 ± 0.05 and k 2 — 1 ± 0.1 gives 
a very good fit to the oblique spindown for all inclina- 
tions. This is roughly consistent with the estimate by 
Gruzinov (2005b). The inferred magnetic field at the 
magnetic equator of a standard neutron star is then 
= 2.6 x 10 19 (PP) 1 / 2 (1 + sin 2 a)- 1 / 2 G which can 
be upto 1.7 times smaller than the estimate from the 
vacuum formula. The fact that oblique rotators can 
lose upto two times more power and evolve faster in 
P — P space than aligned rotators reinforces the predic- 
tion that near the death line there should be an excess 
of pulsars with smaller inclination angles (Contopoulos 
& Spitkovsky 2005). 





Fig. 2. — Oblique pulsar magnctospherc with magnetic inclination a = 60° in the corotating frame: a) Magnetic fieldlines in the fi-fl* 
plane. Color is the magnetic field perpendicular to the plane; b) same as a) but color represents absolute value of the total current |V X B|. 
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Fig. 3. — a) Slices through the 60° magnetosphere. Shown are fieldlines in the horizontal and vertical plane, color on the vertical plane 
is the perpendicular field, on the horizontal plane - toroidal field. Sample 3D flux tube is traced in white, b) Spindown luminosity in units 
of aligned force-free luminosity as a function of inclination. Triangles represent simulation data, and the line is a fit with eq. (3). 



We have developed a numerical method for evolving 
time-dependent force-free MHD equations and have ap- 
plied it to solving a dynamic pulsar magnetosphere. We 
reproduce the closed-open configuration of CKF for the 
aligned rotator, and extend the solution to the gen- 
eral oblique magnetic geometry, calculating the spindown 
rate as a function of inclination. We find that the shape 
of the magnetosphere is controlled by the physics of the 
Y-point. In order for the closed zone to reach the light 
cylinder reconnection must occur near the Y-point in 
the magnetosphere. If this reconnection is for some rea- 
son impeded, the magnetospheric extent may be smaller 
than the light cylinder. Reconnection of the open field- 
lines must also happen when the pulsar spins down and 



the light cylinder recedes. If this reconnection is ineffi- 
cient due to, for instance, intermittent plasmoid ejection 
or inertial loading, the Y-point may lag the light cylin- 
der and lead to braking index less than 3 (Contopoulos 
& Spitkovsky 2005). Y-point dissipation also leads to 
time-dependent phenomena, which when communicated 
to the inner magnetosphere, may leave an imprint on the 
radio emission (e.g. drifting subpulsc phenomena). Such 
behavior can be modeled in our code with different pre- 
scriptions for resistivity, and further research into current 
sheet physics is needed. Availability of self-consistent 3D 
magnetospheric solutions opens the way to the develop- 
ment of quantitative models of phenomena in the mag- 
netospheres of pulsars, magnetars and accretion disks. 
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